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In this work, we propose a self-consistent minimization procedure for functionals in reduced den- 
sity matrix functional theory. We introduce an effective noninteracting system at finite temperature 
which is capable of reproducing the groundstate one-reduced density matrix of an interacting system 
at zero temperature. By introducing the concept of a temperature tensor the minimization with 
respect to the occupation numbers is shown to be greatly improved. 



INTRODUCTION 

Since 1964, after the pioneering work of Hohenberg 
and Kohn [1], density functional theory (DFT) became 
the standard tool for the calculation of groundstate (gs) 
properties of quantum-mechanical systems. There are, 
however, some physical problems which are difficult to 
address in the framework of DFT. These include the de- 
scription of strongly correlated systems, such as the dis- 
sociation of closed shell molecules into open shell frag- 
ments, and the fundamental gap in Mott insulators. Re- 
cently, a promising alternative to DFT was introduced 
which showed success in various fields, ranging from small 
molecules [2-9] to infinite solids [10-12], including the dif- 
ficult cases mentioned above. This method features the 
one-reduced density matrix (1RDM) as central variable 
and is called reduced density matrix functional theory 
(RDMFT). In the theoretical framework of RDMFT, the 
functional form of the kinetic as well as of the exchange 
energy are known exactly in terms of the 1RDM and 
only the correlation part of the two-particle interaction 
energy has to be approximated. However, a minimiza- 
tion of functionals in RDMFT is complicated by the fact 
that at zero temperature there is no noninteracting sys- 
tem reproducing the 1RDM of the interacting system. 
This is in contrast to DFT where the Kohn-Sham sys- 
tem [13] allows for an efficient self-consistent minimiza- 
tion. Therefore, in RDMFT one usually resorts to direct 
minimization routines. 

In the present work, we show that one can indeed con- 
struct a noninteracting system which reproduces a given 
1RDM to arbitrary accuracy, if one considers this system 
to be in grand canoncial equilibrium at finite tempera- 
ture. We therefore effectively model a zero-temperature 
interacting system by a finite-temperature noninteracting 
one. This allows one to construct a self-consistent Kohn- 
Sham minimization scheme for functionals in RDMFT. 

Capitalizing on the freedom of choice for the tempera- 
ture of the Kohn-Sham system, we will furthermore intro- 
duce the concept of a temperature tensor. This concept 
will later on be shown to greatly improve the performance 
of our minimization procedure. 



We will then argue, why the energy value in a numer- 
ical minimization of a RDMFT functional is not a good 
measure of convergence. As alternatives we introduce 
two convergence measures which rely solely on the func- 
tional derivative of the RDMFT functional w.r.t. the 
1RDM. 

Finally, we will investigate the performance of the new 
minimization scheme by applying a common RDMFT 
functional to LiH. It will be shown that the self-consistent 
scheme is very efficient and avoids conceptual difficulties 
prevalent in many other minimization procedures. 



THEORETICAL FOUNDATIONS 

In this work, we will consider systems governed by a 
Hamiltonian H consisting of the kinetic energy T, the 
external one-particle potential V, and the two-particle 
interaction W: 



H =T + V + W. 



(1) 



A quantum-mechanical system is generally described by 
a statistical density operator (SDO) D being a weighted 
sum of projection operators on the Hilbcrt space under 
consideration 



The 1RDM 7(2;, x'), corresponding to a particular SDO 
D, is defined as 



"f(x, x) = tr{£)'i/) + (a;')'i/'(a;)}, 



(3) 



where {'J'(x)} are the common field operators and the 
variable x denotes a combination of spacial coordinate r 
and spin index a (x = (r, a)). An integration over x is 
therefore to be interpreted as an integration over r and a 
summation over a. By construction, j(x, x') is hermitean 
and can therefore be written in spectral representation 



j(x, x') 



ni<j)*{x')(t>i{x). 



(4) 
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The {(f>i(x)} are traditionally called the natural orbitals 
(NO) and the {rii} are the occupation numbers (ON) 
[14]. The conditions that ensure that a given 7(2;, x') is 
ensemble-N-representable, i.e. that it comes from a SDO 
of the form of Eq.(2), are the following [15]. 

< m < 1 (5) 

J2 n *= N ( 6 ) 

i 

{4>i} is a complete orthonormal set (7) 

The set of all ensemble-N-representable lRDMs is given 
by 



{7(21,2/) 

7(x, x') fulfills Eqs.(5), (6), and(7)} (8) 



which is closed and convex. 

Following from the theorems of Hohenberg and Kohn 
[1] , we know that one can formulate a functional theory of 
the 1RDM for the determination of the gs energy. It was 
Gilbert [16] who showed that this theoretical framework 
is also capable of describing systems subject to nonlocal 
external potentials, a task not possible via DFT. We have 
furthermore shown in [17] that this methodology can be 
extended to the case of quantum-mechanical systems in 
grand canonical equilibrium. 

The functional for the energy £[7] of the interacting 
and for the grand potential f2o [7] of a noninteracting sys- 
tem in grand canonical equilibrium are given as 

E[ 1 ]=T[ 1 ] + V ext h] + W[ 1 ] (9) 

no [7] = t[j] + Vextbf] - uNh] - i/PSofr] (10) 

where 



T[l] = 


[ dx' lim ( - 


V 2 \ 

— ) l\x ,x) 


(11) 


VextiH = 


j dxdx'v ext (x 


x')j(x',x) 


(12) 


N[i\ = 


j dxj(x 7 x) 




(13) 


Sob] = - 


- y^(rtj Inn, - 


\- (1 -n l )ln(l -m)). 


(14) 



The functional W[y] for the interaction contribution is 
not known exactly and has to be approximated in prac- 
tice. The noninteracting grand potential can be written 
solely in terms of the one-particle eigenenergies and the 
ONs as 



i 

— (rii In Hi + (1 - rii) ln(l - nj))Y (15) 



In the context of this work, the question of noninter- 
acting (ni)-V-rcprcscntability, i.e. the question which 
lRDMs correspond to a groundstate or equilibrium of 
a noninteracting system, will become important. The 
sets of all zero-temperature ni- V-representable and finite- 
temperature ni-V-representable lRDMs will be denoted 
by and T^, respectively. In the case of zero temper- 
ature a nondegenerate system assumes a pure ground- 
state and the corresponding noninteracting 1RDM will 
be idempotent. Therefore, is on the boundary of T N . 
We have shown in [17] that the gs-lRDM of a Coulomb 
system is in the interior of T N and, therefore, we cannot 
use a noninteracting system at zero temperature to find 
the minimum of an RDMFT functional. In simple terms: 
The 1RDM of interacting particles is never idempotent 
and, hence, it cannot be represented as the 1RDM of a 
Kohn-Sham-type noninteracting system at zero temper- 
ature. At finite temperature, however, for a noninteract- 
ing system with one-particle eigenvalues {^i}, the ONs 
are given by the Fermi-Dirac distribution [18] which can 
easily be inverted: 



(16) 
(17) 



In 



This implies that all lRDMs in the interior of T are in 
r£. Therefore, for every 1RDM in T N there is a 1RDM 
from arbitrarily close to it which allows the utilization 
of a noninteracting system in grand canonical equilibrium 
in a self-consistent minimization scheme. We emphasize 
the term "arbitrarily close" because pinned ONs (i.e. 
or 1) cannot be reproduced by a system at finite temper- 
ature (see Eq. (17)), but every ON arbitrarily close to 
or 1 can. The error introduced by these pinned states 
therefore becomes arbitrarily small. 



SELF-CONSISTENT MINIMIZATION 

The biggest stumbling stone in the numerical mini- 
mization of RDMFT functionals is the incorporation of 
the auxiliary constraints on the ONs and NOs of the 
1RDM. These are particle number conservation n i — 
N, the fermionic constraint < rii < 1, and most impor- 
tantly, the orthonormality constraints of the NOs. Usu- 
ally, the orthonormality of the NOs will be enforced by 
applying an orthonomalization algorithm to the NOs af- 
ter they have been modified, using the information pro- 
vided by the functional derivatives 8E[y]/84>i. These or- 
thonormalization procedures can change several orbitals 
quite significantly which can lead to a slow convergence 
of the minimization routines. 

The main idea of a self-consistent minimization scheme 
is now to approximate the energy surface £^[7] by a sim- 
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pier one whose minimum, incorporating all auxiliary con- 
straints, can be found easily. In our situation, we take the 
information about the derivatives of E[y] at 7 and con- 
struct an effective noninteracting system in grand canon- 
ical equilibrium whose grand potential functional OoM 
has the same functional derivative in 7. The minimum 
of this energy surface is found by a diagonalization of the 
effective Hamiltonian and an occupation of the new ONs 
according to the Fermi-Dirac distribution. The resulting 
cq-lRDM will then serve as the starting point for the 
subsequent iteration. This method automatically incor- 
porates the constraints on the ONs and NOs and we will 
not have to apply subsequent orthonormalizations and 
the like. The success of this scheme, of course, relies on 
the similarity of the energy surfaces of E[j] and fioH- 

We will now proceed to derive the variational equa- 
tions, guiding the determination of 7. 



Effective Hamiltonian 

The effective noninteracting system is constructed such 
that the derivatives of the interacting as well as of the 
noninteracting functional (Eqs. (9) and (10)) coincide. 



«o[7] 



6E[<y] 



Sj(x,x') 5j(x,x') 



(18) 



Because of the possibility of pinned states, this equa- 
tion does not have to be fulfilled exactly. Therefore, 
as mentioned before, our minimization routine may not 
reach the exact minimum but will approach it arbitrarily 
closely. The effective Hamiltonian in spatial representa- 
tion then becomes 

/i e// [ 7 ](x,x') = t[ 1 ]{x,x') + v ext {x 1 x') + 

li5{x - x') + 1/Pa[j\(x, x') + v w [j](x, x'). (19) 

The functional derivatives are given by 

6W[j] 



v w [i\(x,x') 
a[i\{x,x') 



5j(x, x') 

SS [-y] 
Sj(x, x') ' 



(20) 
(21) 



We want to use the chain rule for the functional deriva- 
tive. We therefore need the derivatives of the ONs and 
NOs with respect to 7. They can be obtained using first- 
order perturbation theory, yielding 



5n k 



5~/(x' , x) 
Hk(y) 
5j(x', x) 



Ijtk 



rik - ni 



i(y) 



5~({x',x) ^ rik-ni 



(22) 
(23) 

(24) 



In the following, it will be useful to work in the basis of 
NOs. An arbitrary function g(x,x') is then represented 
by gij, where 



gij = / dxdx'(f>*(x)g(x,x / )(j)j(x'). 



(25) 



The matrix elements /if/ of the effective Hamiltonian 
then become 



V 9m 



1i 

ft 



<t>t(y)), (26) 



where the entropic contribution <Tj is given by 



dS [7] , / 1 - n, 
a t = — = In 



(27) 



The offdiagonal elements are exactly the ones Pernal [19] 
derived in her approach for the derivation of an effective 
potential for RDMFT. They are also simply related to 
the ones Piris and Ugalde [20] introduced in their method 
for an orbital minimization. It has to be noted, however, 
that in our approach the diagonal elements are not free 
to choose but are determined by the thermodynamic en- 
semble. The temperature of the Kohn-Sham system has 
no physical meaning and can be varied to influence the 
convergence behaviour of the minimization routine. If /3 
was small, i.e. if the corresponding effective temperature 
was high, the diagonal part of H e ff will be bigger com- 
pared to the offdiagonal parts. Therefore, after a diago- 
nalization of H e ff , the orbitals will change less. When 
considering the change in ONs, one can investigate the 
limit of /3 — > 0. The diagonal of H e ^ will then just con- 
tain the entropic contribution oi. A solution of Eq. (26) 
will then leave the ONs invariant. We will further inves- 
tigate the behaviour of our self-consistent minimization 
scheme for small j3 later on in this work. 

In the following, we will show how the concept of a 
temperature tensor greatly enhances the adaptability of 
the Kohn-Sham system which will improve the perfor- 
mance of the minimization procedure. 



Temperature tensor 

To understand how the concept of a temperature 
tensor improves the performance of the minimization 
scheme, the following considerations will be helpful. In a 
self-consistent minimization scheme, for a given 1RDM, 
we construct a known (noninteracting) functional whose 
first derivative coincides with the one from the interact- 
ing functional. For a fixed /?, the parameter /x is deter- 
mined by the requirement of particle number conserva- 
tion. j3 can then be varied to modify how narrow the 
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(a)First ON, ai 



50 



0.2 0.4 0.6 O.i 

(b)Second ON, a 2 = 1 



FIG. 1: Projected energy surfaces for the model of Eq. (28) 
with qi = 50, a 2 = l,/3 = 0.11. 



nonintcracting energy surface should be. However, sec- 
ond derivatives with respect to the ONs may differ quite 
substantially and a value of (3 which describes the en- 
ergy surface w.r.t. one ON well might describe others 
quite badly. A simple example is the following quadratic 
two-state model functional E[ni,n 2 ] without orbital de- 
pendence. 

E[m, n 2 ] = y(m - 0.5) 2 + ^ (n 2 - 0.5) 2 (28) 



Ei[n{\ + E 2 [n 2 }. 



(29) 



eff 



The choice of ai = 50 and a 2 = 1 leads to h\\ 
ei = -0.225 + /i and h% s = e 2 = 0.00450 + fx in Eq. 
(26). The corresponding projected grand potential sur- 
faces from Eq. (15) are plotted in Figure 1 for (3 = 0.11. 
As one can see, the choice of P = 0.11 models the first 
energy surface quite well, but the second one fails to be 
reproduced. One would like to have some sort of state- 
dependent (3i which can be related to the second deriva- 
tives. However, before one can use such a construct, one 
has to confirm that it corresponds to an energy-surface 
whose minimum can easily be found. 

We arc now going to show that this is possible by a 
slight variation of the definition of grand canonical en- 
sembles. We consider the following generalized SDO- 
grand potential functional 



&[D] = tr{D{B{H - fxN) + hit))}, 



(30) 



where B is an arbitrary hcrmitean operator on the Fock- 
space. The same proof as in [21] now leads to the follow- 
ing variational principle 

G[D]> <6[D eq ], (31) 

where the equality is only fulfilled if D = D eq , with 



D eq = e-^-^/Z e 



tr{e 



-B(H-p.N) 



}• 



(32) 
(33) 



For a nonintcracting Hamiltonian and a B for which 




(a)First ON, /3 = .5 



(b)Second ON, 2 = 5.55 



FIG. 2: Projected energy surfaces for the model of Eq. (28) 
with Qi = 50, a?2 = 1- The choice of r\ — 0.5 leads to /?i = 0.11 
and @2 = 5.55. 



[£>, H] = the Fermi Dirac relation reads 

1 



e ft( £i -M) + i 

1 ( 1 - n, 
Si - [i = — In 



Pi 



n: 



(34) 
(35) 



where Pi denotes the i-th eigenvalue of B. This leads to 
the following expression for the grand potential 



n 



i 

— (ni Inn, + (1 - rii) ln(l - rii))) (36) 

Pi ' 



= y^Ooi[nj, Pi 



(37) 



Where in the case of a scalar temperature we just had 
one parameter to construct our effective noninteracting 
system, we now have one for each ON. A straightformard 
utilization of this freedom would be to let the second 
derivatives of the energy functional with respect to the 
ONs of the interacting functional and the noninteracting 
one be proportional to each other. 



dnf 
1 



71,(1 - Hi) 



d 2 E[ 1 ] 



(38) 
(39) 



where 77, the proportionality factor, is the only global 
parameter. In our model (Eq. (28)), this yields 



P, 



V 



1 



ai n,(l 



(40) 



i] = 1 lets the second derivatives of interacting and non- 
interacting functional be equal whereas an increase (de- 
crease) of 77 leads to a spreading (compression) of the 
noninteracting energy surface. As can be seen from Fig- 
ure 2, with a good choice of 77 (in our model r\ = 0.5) 
one can reproduce the different energy surfaces simulta- 
neously. 



5 



1.) diagonalize ^ k \x 








2.) from [ni k) \ and 


(j>[ k) \ calculate H e ^^ by Eq. (26) 






f 

3.) diagonalize # e -W*) to get 


a new set of and 






j 


* 


4.) construct |^ fe+1 ^ 


> using Eq. (35) and find such 


that ™, <t+1) = N 






I 


5.) from {<if' +1) (s)} 


and (nf +1) ) build new 1RDM 








1 


6.) mixing of y(x,%') with yields *y( k+1 ) (x,x') 



FIG. 3: Self-consistent minimization scheme in FT-RDMFT 



We can now construct a self-consistent scheme for the 
minimization of E[j] which we sketch in Figure 3. A 
mixing of lRDMs is straightforward, because T N is a 
convex set. 



Small step investigation 

We showed in the previous considerations that one can 
employ the Kohn-Sham system in FT-RDMFT to con- 
struct a self-consistent minimization scheme. However, 
this does not ensure that an application of this scheme 
will actually lead to a minimum of the functional. This is 
a common problem of minimization schemes, but in the 
following we are going to show that for small steplengths 
our method will definitely lead to a decrease of the value 
of the functional under consideration. As we argued be- 
fore, choosing a smaller j3 will lead to smaller changes 
in ONs and NOs. Starting from a given 1RDM 7, we 
therefore apply first-order perturbation theory to get the 
modified 1RDM 7'. By the virtue of Eq. (26), 7 leads to 
the effective Hamiltonian H e ^ . A diagonalization under 
the assumptions of first-order perturbation theory then 
yields the following new eigenvalues e\ and eigenstates 



The new ONs, resulting from our modified cigcncncrgies, 
become 



1 



1 _)_ e Pi(e'i-V-A-n) ' 



(43) 



where one had to introduce the chemical potential cor- 
rection A M to ensure particle number conservation. With 
Eqs. (26) and (27) one gets 



(44) 



rij + (1 - rii)e 



Pi 



:-A„ 



Expanding Eq. (44) in orders of fa and retaining the 
leading contribution, we get 



Srii = n. 



(45) 
(46) 



This result is very similar to the steepest-descent method 
with an additional factor of n, (n; — 1). This additional 
term tries to keep the ONs in the allowed set < n» < 1. 
A^ ' can now be found by the requirement of particle 
number conservation, 



A(0) = E i A^ z i]gg 



(47) 



The overall change in the 1RDM up to first order in fa 
is then given by 



A 7y = in - Sijfii 
The energy changes accordingly as 



AE 



dxdx 



-A<y(x',x) 



5j(x, x') 
^— ' On, ^— ' <=■-■ — «•-• l -> 



(48) 
(49) 

(50) 
(51) 

(52) 



AE! 



AE 2 



We see that the energy change AE seperates into two 
parts. One is determined by the change in ONs, and 
one comes from the change in NOs. In the following we 
are going to investigate these two different contributions 
seperately. 



eff 



(41) 

Occupation number contribution 



<j>i{x) = (f>i{x) + ^2 J ^ £ 0j( x ) ( 42 ) We will now show that the first term in Eq. (52), which 

1 3 is due to the change in ONs, is negative for appropriately 
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small stcplcngths. 



(53) 



(54) 



For brevity, we introduce Cj = ^'^g. „"('„. ~T) ■ Using Eq. 
(47) then leads to 



A£ 1 = ^^n 3 (n,-l) 



E c > 



• (56) 



Because every ON n, fulfills < < 1 and every is 
greater 0, this leads to the conclusion 



A£a < 0. 



Natural orbital contribution 



(57) 



We can now turn to the second term in Eq. (52) which 
represents the energy change due to the change in NOs. 



By using Eq. (35) this transforms to 



A£, 



n e //|2 



'^J /3 



i i n f _ i i n f y 



(58) 



(59) 



We see that for an arbitrary choice of we cannot en- 
sure the negativity of Af2 2 - But if we use a constant /3, 
we get 



which is nonpositive for < < 1: 
A£ 2 < 0. 



(61) 



We have shown that for small enough the ON 
change will always decrease the grand potential, regard- 
less of wheather one chooses a constant temperature or 



a temperature tensor. When considering changes in the 
NOs, one has to fall back to constant temperature to en- 
sure a decrease of the functional value. We will use these 
findings in our numerical implementation of the mini- 
mization scheme as demonstrated later on in this work. 



Convergence measures 

We have now all the neccessary tools at hand to iter- 
atively minimize a functional £[7]. We need, however, 
some measures to judge if a calculation is converged. 
There are two main reasons why using the energy it- 
self as convergence measure is disadvantageous. Firstly, 
often the calculation of derivatives is not accurate and 
consequently a derivative-based minimization may lead 
to a fixpoint where £[7] is not minimal. This leads to a 
sign change of the convergence measure and implement- 
ing the strict decrease of energy as a requirement of the 
minimization procedure will then lead to a starting point 
dependent result. Secondly, because the true minimal en- 
ergy is not known, one would have to judge convergence 
from the change in £'[7] after iterating the minimization 
routine, i.e. a small change in £[7] indicates a relative 
closeness to the real minimum. This might pose a prob- 
lem if the minimum of the energy-surface, as defined by 
the derivatives of £[7], is very shallow or, worse, if the 
minimization procedure leads to a slow approach to the 
minimum. An example for such a situation is discussed 
in the summary. 

Because of these problems, we would rather use a 
strictly positive convergence measure which goes to 
if the 1RDM approaches the minimum of the energy- 
surface, as defined by the derivatives of £[7]. We will 
establish our choice of convergence measures on the fol- 
lowing two observations. 

• In the minimum, the derivatives with respect to the 
ONs will be equal for unpinned states. 

• In the minimum, the effective Hamiltonian H e ^ 
will be diagonal. 

The first observation allows us to define a convergence 
measure \n f° r a minimization with respect to the ONs: 



1 



N, 



unpinned 



N unpinri 

E 



dE 
dm 



N, 



unpinned 



PT cd dE_ 
dm 



E 



(62) 



(63) 



The second statement leads to the following definition 
of x% as a convergence measure for a minimization with 
respect to the NOs: 



1 E 



N \h e ff\ 2 



N - 1 



2 

i w i 



(64) 
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If a minimization is converging, both measures should 
approach 0. 

EXAMPLE 

We test the self-consistent procedure for the case of 
solid LiH at zero temperature by using the FP-LAPW 
code Elk [22]. The exchange-correlation energy will be 
modelled by the a functional E" c [y], as introduced in 
[12]. The interaction functional W[j] is then given as a 
sum of the classical Hartree energy functional Eh [7] and 
the exchange-correlation functional E" c [7] . 

*E c [7] = ~I>? n i / dxdx ' 

w(x, x')fa* (x')fa(x)<p* (x)fa (x) (65) 

We choose this functional because it exhibits several 
properties making it difficult to be minimized. It will 
lead to several fully occupied, i.e. pinned states. There- 
fore, as argued before, there is no noninteracting system 
at finite temperature reproducing this 1RDM exactly but 
there will be one leading to an equilibrium 1RDM ar- 
bitrarily close. The minimization of the a functional is 
therefore a good test for the minimization scheme leading 
to boundary minima on T N . Furthermore, the a func- 
tional exhibits divergencies in the derivatives w.r.t. the 
ONs for rii — > 0. If, in the minimum, there are ONs close 
to (and there will be if one considers enough NOs) this 
might lead to convergence problems of the minimization. 

We will now investigate the performance of the self- 
consistent minimization scheme w.r.t. ON- and NO- 
convcrgcncc. It turned out in the course of our inves- 
tigations that a mixing of 1RDMS, according to point 6 
in Figure 3, does not improve our results and we therefore 
abstain from it. 

Occupation number minimization 

We have minimized the a functional for a — 0.565 with 
three methods. First, we have used the steepest-descent 
method, as implemented in Elk. The second method is 
the self-consistent FT-RDMFT minimization with con- 
stant /3, and finally we have employed a temperature 
tensor j3i of the form of Eq. (39) with parameter ?/. In 
all three methods, we chose all parameters to achieve 
fastest convergence. The results, which are shown in 
Figures 5(a)-5(f), show that both self-consistent Kohn- 
Sham minimizations lead to a faster convergence than 
steepest-descent. A dramatic improvement in the speed 
of convergence is achieved by employing a temperature 
tensor. The slow decrease of Xn m Figures 5(d) and 5(f) 
for the steepest-descent and constant-/? methods can be 
attributed to the following fact. For these two methods, 



One full effective Hamiltonian diagonalization 
with constant ft 



ON minimization with /Jj-tensor 



FIG. 4: Full minimization scheme 

the ONs, which will be pinned at the equilibrium, ap- 
proach their final values quite slowly. Therefore, their 
derivatives contribute to x\ YY& Eq- (62) even after sev- 
eral iterations. 

Full minimization 

We can now turn to the problem of minimizing E[y] 
with respect to both ONs and NOs. We find that the 
overall performance of this full minimization is greatly 
improved by introducing a ON-minimization after ev- 
ery NO-minimization step (see Figure 3). Because we 
have seen in the previous section that this can be done 
very efficiently, this increases the runtime of a full mini- 
mization run only negligibly. The deeper reason for the 
improvement of the convergence by inclusion of an ON- 
minimization is the following: It typically happens that 
two states fa and faj have similar eigenvalues in but 
considerably different ONs. A diagonalization of H e ^ 
then yields a strong mixing between these states. If the 
ONs were not updated, one might be led away from the 
minimum of the grand potential functional. A subse- 
quent ON-minimization remedies this problem and as- 
signs the optimal ON for each NO. We show a sketch of 
the full minimization scheme in Figure 4. An applica- 
tion of this scheme to LiH then leads to the results de- 
picted in Figure 6. Again, we see a tremendous increase 
in speed and accuracy for the self-consistent Kohn-Sham 
minimization scheme compared to the steepest-descent 
method. The steepest-descent method shows a very slow 
convergence, which can be attributed to the orthonor- 
malization of NOs. The increase of the energy curves in 
Figures 6(c) and 6(e) is due to the approximative nature 
of the derivatives. 
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FIG. 5: Energy E and ON-convergence measure Xn for ON minimizations of the a functional, with a — 0.565, applied to LiH. 
The red, short dashed lines stand for a steepest-descent minimization, the blue, solid ones for a sc-Kohn-Sham minimization 
with constant /?, and the black, long dashed ones for a sc-Kohn-Sham minimization with adaptive /3i . r denotes the parameter 
value for taurdmn in Elk, whereas /3 and rj are defined via Eqs. (10) and (38). 
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FIG. 6: Energy E and NO-convergence measure x<P f° r NO minimizations of the a functional, with a = 0.565 applied to 
LiH. Both variables are plotted against the number of NO changes. After each change in NO there follows a complete ON 
minimization. The red, dashed lines stand for a steepest-descent minimization whereas the blue, solid ones depict a sc-Kohn- 
Sham minimization with constant /3. r denotes the parameter value for taurdmc in Elk, whereas /? is defined via Eq. (10). 
The increase of energy is due to the fact that the derivatives are calculated only approximately. 
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FIG. 7: Energy convergence for different effective tempera- 
tures. 



FIG. 8: NO-convergence for different effective temperatures. 



SUMMARY AND OUTLOOK 

In the present work, we have introduced a self- 
consistent Kohn-Sham minimization scheme in the theo- 
retical framework of RDMFT. We have defined measures 
which allow us to judge the convergence of a calculation 
without having to resort to the energy. We could show 
that this self-consistent procedure is superior in many 
respects compared to the steepest-descent method, espe- 
cially considering a minimization w.r.t. the NOs. The 
important parameter in the minimization scheme is the 
effective temperature (3 and the speed of convergence cru- 
cially depends on it. In Figures 7 and 8 we show the 
behaviour of the minimization scheme for three different 
choices of (3. 13 = 1 represents the optimal value, i.e. the 
value for which the convergence measure x| decreases the 
fastest. We see that the energy reaches its fixpoint after 
approximately 300 iterations. An increase of /3 to /? = 2 
seemingly speeds up the energy convergence, but from 
one can see that after about 100 iterations the minimiza- 
tion fails to diagonalize h e ^ any further. The changes in 
the 1RDM, whose amplitudes are determined by /3, be- 
come too big and the 1RDM jumps around the fixpoint 
of the energy. Without considering this would have 
been difficult to detect which illustrates the importance 
of a convergence measure which is independent of the en- 
ergy value. One might argue that this choice of (3 still 
leads to a fixpoint very close to the optimal one, but this 
cannot be ensured for all problems and all choices of (3 
and therefore has to be seen in the actual example as 
rather accidental, i.e. fortunate. A further increase of f3 
to f3 — 5 then exposes this problem more dramatically. 
The energy apparently reaches a fixpoint. But this fix- 
point is considerably above the optimal one. Just having 
the energy at hand, this would have been difficult to de- 
tect. But directly shows that the minimization is far 
from being converged. 

One important feature, which can be extracted from 



Figures 7 and 8, is that all three parameters lead to 
a similar energy vs. iteration curve. Apparently, a 
minimization-run with (3 being too big is able to lead 
to the vicinity of the fixpoint. An utilization of this fact 
would now be to use an adaptive (3 rather than a con- 
stant one. One could start with a big /? till the energy 
does not change anymore and than decrease (3 until 
surpasses the convergence threshold. 

We expect that our successful demonstration of an ef- 
ficient minimization scheme in RDMFT will support the 
investigation and development of functionals and there- 
fore encourage further work in this field of research. 
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